rm(list=ls())

packs <- c("tidyverse", "ggplot2", "stargazer", "sp", "maptools", "rgeos", "rgdal") 
# install.packages(packs, dependencies = TRUE) ## uncomment if you need to install any of the necessary packages
lapply(packs, library, character.only = TRUE)

setwd("C:/Users/mar20012/Dropbox/AMR_1948/Depopulation_Social_Cohesion/IS_Supplementary_Materials")
#setwd("~/Dropbox/AMR_1948/Depopulation_Social_Cohesion")


village_master <- read_csv("data/village_sample_initial.csv")

df <- readRDS("data/df.rds") %>% 
  as_tibble()

################################################################
################################################################


################################################################
################################################################
## Figure 1: Map of Village Evacuation Outcomes
################################################################
################################################################

#crs is 3857. lat/long need to be imported at 4326

villages <- readOGR("data/shapefiles/villages_id_crs3857.shp")
un_border <- readOGR("data/shapefiles/un_border_crs3857.shp")
roads <- readOGR("data/shapefiles/roads_crs3857.shp")
outside_border <- readOGR("data/shapefiles/outside_border.shp")

# prepare the 3 components: coordinates, data, and proj4string
coords <- village_master[!is.na(village_master$latitude) & !is.na(village_master$longitude) & village_master$refugee_camp == 0 & village_master$bedouin == 0& village_master$pop0 == 0 & village_master$jewish_only == 0, c("longitude", "latitude")]

# filter out villages with no population (pop0 = 1), refugee = 1, bedouin = 1
data   <- dplyr::select(filter(village_master, !is.na(latitude), !is.na(longitude), refugee_camp == 0, bedouin == 0, pop0 == 0, jewish_only == 0),-latitude, -longitude)
crs    <- CRS("+init=epsg:4326")


# make the spatial points data frame object
village_master_shape <- SpatialPointsDataFrame(coords = coords,
                                               data = data, 
                                               proj4string = crs)
village_master_shape <- spTransform(village_master_shape,CRS("+init=epsg:3857"))


village_master_shape_sf <- sf::st_as_sf(village_master_shape)


ggplot() + 
  geom_polygon(data = outside_border, aes(x = long, y = lat),colour = "black",fill=NA)+
  theme_void()+
  geom_path(data = un_border, aes(x = long, y = lat,group=id),colour = "gray20",fill=NA,alpha=.5)+
  geom_path(data = roads, aes(x = long, y = lat,group=id),colour = "gray40",fill=NA,alpha=.25, lty=2)+
  geom_sf(data=village_master_shape_sf, aes(shape = factor(depop_outcome), color = factor(included)),size=1) +
  #geom_point(data=village_master_shape, aes(x = longitude, y = latitude, shape = factor(depop_outcome), colour = factor(included)),size=1) +
  scale_shape_discrete(name="evacuation outcome")+
  scale_colour_manual(name="areas of operation", values = c("black","grey50")) +
  theme(legend.position=c(.2, .8),
        legend.title = element_text(colour="black", size=12, face="bold"),
        legend.text = element_text(colour="black", size = 10, face = "bold"))#+

ggsave(paste("figures/Figure1_Evacuation_Outcome_Map.pdf", sep = ""), width=4.83, height = 10, units = "in")


################################################################
################################################################
## Figure 3: Distribution of Village Evacuation Outcomes
################################################################
################################################################

df_fig3 <- data.frame(main_var = c("preemptive evacuation", "remain", "remain")
                  , factor_var = c("preemptive evacuation", "reactive evacuation", "no evacuation")
                  , count=c(nrow(dplyr::filter(df, depop_preemptive == 1))
                                 , nrow(dplyr::filter(df, depop_violent == 1))
                                 , nrow(dplyr::filter(df, depop_preemptive == 0 & depop_violent == 0))
                            )
                  , order = c(3, 2, 1)
                  , order2 = c(3, 2, 1)
)

df_fig3 <- data.frame(main_var = c("preemptive evacuation", "remain", "remain")
                      , factor_var = c("preemptive evacuation", "reactive evacuation", "no evacuation")
                      , count=c(nrow(dplyr::filter(village_master, v_survey_collected == 1, depop_outcome == "preemptive evacuation", included == "contested"))
                                , nrow(dplyr::filter(village_master, v_survey_collected == 1, depop_outcome == "reactive evacuation", included == "contested"))
                                , nrow(dplyr::filter(village_master, v_survey_collected == 1, depop_outcome == "no evacuation", included == "contested"))
                      )
                      , order = c(3, 2, 1)
                      , order2 = c(3, 2, 1)
)


ggplot(data = df_fig3, 
       aes(x = reorder(main_var, order), y = count, fill = reorder(factor_var, order2))) + 
  geom_bar(stat="identity", position="stack") +
  labs(x = "evacuation outcome" , y = "number of villages") +
  scale_fill_grey() + 
  guides(fill=guide_legend(title="evacuation outcome")) +
  theme_light() +
  theme(legend.position=c(.775, .825), panel.grid.major = element_blank(), panel.grid.minor = element_blank())





ggsave(paste("IS_Supplementary_Materials/figures/Figure3_Evacuation_Freq.pdf", sep = ""), width=4.83, height=4.75, units = "in", pointsize=11)


################################################################
################################################################
## Figure 4: Distribution of the Social Cohesion Index
################################################################
################################################################

ggplot(data=df,aes(social_cohesion_add))+
  geom_bar()+
  xlab("social cohesion index") +
  theme_classic()+
  ylab("number of villages") + 
  #ggtitle("Distribution of Social Cohesion Index Values") +
  theme_light() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

ggsave(paste("IS_Supplementary_Materials/figures/Figure4_SocialCohesion_Freq.pdf", sep = ""), width=4.83, height=4.75, units = "in", pointsize=11)


################################################################
################################################################
## Figure 5: Map of Evacuation Outcome x Social Cohesion
################################################################
################################################################
village_master_shape_plus_soc <- readRDS("paper_input/village_master_shape_plus_soc.rds")

ggplot() + 
  geom_polygon(data = outside_border, aes(x = long, y = lat),colour = "black",fill=NA)+
  theme_void()+
  geom_path(data = un_border, aes(x = long, y = lat,group=id),colour = "gray20",fill=NA,alpha=.5)+
  geom_path(data = roads, aes(x = long, y = lat,group=id),colour = "gray40",fill=NA,alpha=.25, lty=2)+
  geom_point(data=village_master_shape_plus_soc, aes(x = longitude, y = latitude, color = social_cohesion_add, shape=factor(depop_outcome,labels=c("preemptive evacuation","no evacuation","reactive evacuation"))))+
  scale_shape_discrete(name="evacuation outcome")+
  scale_colour_gradient(name="social cohesion index", low = "gray", high = "black")+
  #scale_colour_discrete(name="social cohesion index", c("gray25", "gray50", "gray75", "black"))+
  theme(legend.position = c(0.2, .8)) #+ 
  #ggtitle("Community Evacuation and Social Cohesion")

ggsave(paste("IS_Supplementary_Materials/figures/Figure5_Evacuation_and_SocialCohesion_Map.pdf", sep = ""), width=4.83, height=10, units = "in", pointsize=11)

################################################################
################################################################
## Figure 6: Distribution of Village Evacuation Outcomes x Social Cohesion
################################################################
################################################################


df_fig6 <- data.frame(main_var = c("0", 
                               "0",
                               "1",
                               "1",
                               "2",
                               "2",
                               "3",
                               "3")
                  ,  factor_var= c("remain",
                                   "preemptive",
                                   "remain",
                                   "preemptive",
                                   "remain",
                                   "preemptive",
                                   "remain",
                                   "preemptive")
                  , count=c(nrow(dplyr::filter(df, social_cohesion_add == 0, depop_preemptive == 0,district!="Tulkarm"))
                            , nrow(dplyr::filter(df, social_cohesion_add == 0, depop_preemptive == 1,district!="Tulkarm"))
                            , nrow(dplyr::filter(df, social_cohesion_add == 1, depop_preemptive == 0,district!="Tulkarm"))
                            , nrow(dplyr::filter(df, social_cohesion_add == 1, depop_preemptive == 1,district!="Tulkarm"))
                            , nrow(dplyr::filter(df, social_cohesion_add == 2, depop_preemptive == 0,district!="Tulkarm"))
                            , nrow(dplyr::filter(df, social_cohesion_add == 2, depop_preemptive == 1,district!="Tulkarm"))
                            , nrow(dplyr::filter(df, social_cohesion_add == 3, depop_preemptive == 0,district!="Tulkarm"))
                            , nrow(dplyr::filter(df, social_cohesion_add == 3, depop_preemptive == 1,district!="Tulkarm"))
                            
                  )
)


ggplot(data = df_fig6, 
       aes(x = main_var, y = count, fill = factor_var)) + 
  geom_bar(stat="identity", position=position_dodge()) +
  labs(x = "social cohesion index" , y = "number of villages") +
  scale_fill_grey(labels = c("preemptive evacuation", "remain")) + 
  guides(fill=guide_legend(title="evacuation outcome")) +
  theme_light() +
  theme(legend.position=c(.775, .825), panel.grid.major = element_blank(), panel.grid.minor = element_blank())


ggsave(paste("IS_Supplementary_Materials/figures/Figure6_Evacuation_by_Cohesion.pdf", sep = ""), width=4.83, height=4.75, units = "in", pointsize=11)
